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1. Introduction 



O I The Jacobi Last Multiplier (JLM) [1][7] plays, in first order linear partial differential 

^ I equations, a role similar to the integrating factor in first order ordinary differential 

equations. If one can guess a JLM, it is possible to find the general solution of the 
_ equation, or for equations with more than two variables to reduce the number of 

I variables. However in the case of quasilinear first order partial differential equations 

^ ■ we can always integrate them going over to the characteristics. So the role of JLM 

is certainly not crucial as an integrating tool of PDEs, but it has an important role 
in many other applications. For example, JLM has recently received a great deal of 
attention in the theory of A-symmetries of differential equations [9l[15] . 

In this work, we present a first approach to an equivalent concept for difference 
equations. In the case of partial linear difference equations no integration technique 
equivalent to the use of the characteristics exists pill] so the use of the JLM can be very 
proficuous. 

Section 2 is devoted to a short review of JLM in first order partial differential 
equations, in particular the method to obtain the equation satisfied by a JLM. In section 
3, the case of a difference equation in a two dimensional lattice is fully developed together 
with some examples. The conclusions presented in Section 4 are devoted to a summary 
of the results obtained, showing the difficulties which appear when extending the method 
to a higher number of variables, and some future perspectives. 
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2. The continuous Jacobi last multiplier: a review 

Let us consider a first order linear partial differential equation: 

N 

xu = o, x = 5^/wa,„, 



(1) 



i=l 



where f^^^ are some smooth functions on the variables x^^\ If — 1 functionally 
independent solutions of this equation are known, u^-^\ . . . ,u^^~^\ we can write the 
following determinant 

/ 9m du \ 



X 



det 



V 



duW 



(2) 



for any function u{x^^\ . . . , x^"^). The determinant is zero only if m is a solution of 
equation ([1]), as it can be easily shown since the elements of the matrix are the coefficients 
of the forms du^'^^ and the functions u, u^^\ . . . , u^^~^^ are functionally dependent if u is 
a solution of the differential equation. Then the equations 

have the same set of solutions and must be proportional as — 1 functionally 
independent solutions fix the coefficients of the linear first order PDE up to a global 
factor: 



The multiplicative factor M is called the Jacobi last multiplier. 

If we expand the determinant ([2]) using the first row, we find that equation 
be written as 

N 



du 



du 

dx("-' 



i=l 



du 



(4) 



can 



(5) 



where A^'^^ are the corresponding minors of the first row of the matrix in ([2]). Comparing 
the coefficients of the derivatives of in ([S]), we get 

= M/('=), k = l,...,N. (6) 

We can write now a differential equation satisfied by the functions A^''^ 

/ 9^(1) duW duW duW \ 

I dxW ' ' ' axC^-i) ■ ■ ■ a^^cn) \ 

(-l)'"Met 



V 



dxW 

a„(jv-i) 



A'^''\ k 



, N. (7) 



This is a general relation arising from the particular form of the functions A^''\ written 
as minors of the matrix ([2]), and it does not depend on the fact that the functions m^^-* 
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are solutions of the linear homogeneous first order PDE ([T]) with coefficients f^''\ In 
fact, the equation for A^^^ can be considered as a consistency condition for the system 
of equations in u^^^ with k = 1, . . . , N, j = 1, . . . , N — 1. 
Let us write ([7]) as 

= det(Bi, ...,Bk,---,BN), ^ (8) 

where by Bi we denote the z-column of the matrix ([2]) and the symbol means that 
the corresponding /cth-column of ([2]) is removed. Deriving ([8]) with respect to x^''^ and 
summing over k we get: 

^ ^ Fib ^ 

E7M^= E (-irMet(5,,...,— (9) 



Since 



-1 i,k=l,ij^k 



dB, f 92^(1) a^n^^-^) \^ dBk 



(10) 



we have 

det(i?i, . . . , -Q^^ • • • , -Sfc, • • ■ , -Bat) = det(i?i, . . . , • • • i -Bfc, • • • , ^at) 

= (-1)^=— 1 det(Si, . . . , i?;v). (11) 

Finally, since (— 1)*^(— 1)'^~*~^ = —(—1)* the sum of all terms in ffTOj) is zero and 
consequently we get the equation: 

+ ^^^^ 
In three dimensions, this expression is the classical formula of vector calculus stating 

that the divergence of the cross product of two gradient vectors is zero. In an arbitrary 

dimension it can be written using exterior products and differential forms (see for 

instance [B]). 

We can derive an equation for M, differentiating ([6]) with respect to x^^"^ and 
summing over all k between 1 and A^. We get 

^ = ^X^afw ^^•^^''^a^ ^^^^ 
fc=i fc=i fc=i 



and consequently, taking into account (jT2l) . we obtain: 

fe=i fc=i 

This equation depends only on the differential equation ([T]) and thus M does not 
depend on any particular solution. From a practical point of view, equation ( IT4|) is 
an inhomogeneous version of the original equation ([T]). However, as in the method of 
integrating factors, if we know a particular solution of f lT^ . we could use it to compute 
a solution of ([1]). This is exploited in the examples presented in the next subsection 
where we consider a few examples and we compare the results obtained by the use of 
JLM with those obtained through other methods of integration of the linear first order 
partial differential equations. 
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2.1. Examples 

As a simple illustration of the method, let us consider the following examples of partial 
differential equations in two and three independent variables. 

2.1.1. Two independent variables 

yu^ + xuy = 0. (15) 

From ([T3j) the equation for M is: 

ydJogM + xdylogM = 0. (16) 

An obvious solution of this equation is M = 1. Then, 

= M/i = y, = Mf2 = x (17) 

and ([6D reduces to the compatible overdetermined system of equations: 

= -X, Uy = y. (18) 

Solving this system we find a non trivial solution of the original partial differential 
equation 

u{x,y) = ^{y''-x''). (19) 

The general solution can be obtained by computing the characteristic variable ^ = y'^—x'^ 
and is given by: 

u{x,y) = Fiy^-x''), (20) 

where F is an arbitrary function of its argument defined by the initial conditions or 
from the boundary values. 

2.1.2. Three independent variables 

x{x + y)u^ - y{x + y)uy + z{x - y)u^ = 0. (21) 

The equation satisfied by a Jacobi last multiplier M is 

d d d 

x{x + y) — log M - y{x + y)—\ogM + z{x - y)—\ogM + 2{x - y) = 0. (22) 

Looking for a particular solution of this equation, for instance M = M{z), we find 
M = ^, and the system of equations we have to solve is (with u^^^ = u, u^'^^ = v): 

x{x + y) y{x + y) x-y 

z z z 

Given a solution u of (12T|) . fl23|) is an overdetermined system for v. For instance, we can 
take u{x, y, z) = xy and the new solution v is obtained from the following overdetermined 
system of equations 

X + y X — y 

Vz = —, yVy - XV^ = . (24) 
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A solution of (jM]) is 

. = (25) 
z 

From the method of characteristics we find that any solution of fl2T|) is a function of the 
two particular solutions we have found 

u{x,y,z) = F (^xy,^^—^^ . (26) 
3. Difference equations 

A difference equations for one dependent variable m is is a relation between the 
function in various points of a lattice. If the lattice is r dimensional it can be put 
in correspondence with the points of an r-dimensional space. 

An ordinary difference equation (OAE) is a difference equation on a one dimensional 
lattice. In this case the lattice is given by an ordered sequence of points on a line 
characterized by their relative distance (see Fig.l). If Xi and Xj+i are two subsequent 
points, their distance — Xj| will be hi. We can then introduce a x-shift operator 
Tx such that T^Xi = Xj+i and in term of it we can construct delta operators which in 
the continuous limit, when /ij — )■ 0, go over to the derivative. An example of such delta 
operator is given by the right shifted discrete derivative 

, , u(xi+i) — u(xi) (Tx — l)u(xi) 
Axu{xi) = ^ +^ ^ = ^ — ^ ^ ' . (27) 

Xi^\ '^i i^i 



In some instances it may be more convenient to introduce symmetric delta operators [13] 
as, for example, 

Xi+l ~ Xi-i lli + lli-l 

An OAE of order n, i.e. involving n + 1 points of the lattice, can thus be written as 

£{xi, Ui, TxUi, T^Ui, ■ ■ ■ T^Ui) = 0, (29) 

or, equivalently, in term of the operator delta as 

T{xi, Ui, AxUi, Alui, ■ ■ ■ A^Mj) = 0. (30) 

However the equation fl29l) (or fl30l) ) is not completely defined unless we specify the 
values of the distance between the various lattice points hj, j = i, ...,i + n involved in 
the equation. This implies that an OAE will be defined only if we attach to it a second 
equation which defines the lattice. The set of these two equations is called a Difference 
scheme. 

In many instances, when the equation comes from some physical problem, the 
lattice is a priori given, as for example, when all the points are equidistant so that 
hi = h. In this case the lattice equation is trivial Xj+i — Xi = h. However there may be 
situations, as, for example, discretizing a continuous differential equation to solve it on 
the computer, when we want to take advantage of the freedom of the lattice to preserve 
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(a) the index line (b) the x-lme 

Figure 1: 1 -dimensional lattice grids 
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(a) the index plane (b) the x, y~plane 

Figure 2: 2-dimensional lattice grids 

in the discretization other properties of the continuous system like its symmetries |14j . 
In such a situation the lattice may be defined by a non trivial equation maybe also 
depending on the dependent variable so as to have a denser grid when the solution 
varies rapidly. 

In the case of OAE's there is at least one natural parametrization of the differences 
which in the continuos limit go to the corresponding derivatives and which simplifies the 
discretization procedure [12]. Such parametrization gives a one-to-one transformation 
between the lattice differences, discrete approximations of the derivatives, and the lattice 
points. 

A similar situation exists also in the case of partial difference equations (PAE's), 
however in this case the definition of the lattice must be given by compatible equations 
as the independent variables depend on several indices (see Fig. 2 for the 2-dimensional 
case where Xn,m and y„,m, depend on two indices). In general the difference scheme will 
be given apart from the P AE, by a set of equations which define the lattice and depend 
on the number of independent variables and on the problem we are solving [13] . However 
in this case, up to now, no natural parametrization exists which in the continuous limit 
goes to the corresponding derivatives and which simplifies the discretization procedure. 
Work on this is in progress [T0l[T6] . 
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The solution of linear OAE's follows the standard technique of solving ordinary 
differential equations. The solution is given by a linear combination with arbitrary 
coefficients of powers of the independent variable and the exponents are defined by a 
characteristic polynomial. In the case of linear PAE's the situation is more complicate 
(as is also the case for partial differential equations). As one can read in [3] The method 
of trial and error is still one of the basic methods for obtaining explicit solutions. If the 
PAE has constant coefficients then a few techniques can be found in Jordan book [8], 
such as Laplace method of generating functions, or the method of Fourier, Lagrange and 
Ellis [5]. If the PAE does not depend explicitly on one of the two independent variables 
then Boole symbolic method can be applied [2]. 

Consequently, it seems particularly important to extend the last Jacobi multiplier 
technique to the case of linear PAE's as this will provide solutions also in the case of non 
constant coefficient PAE's. In the following, for the sake of simplicity, we will consider 
the case of an orthogonal lattice, where the independent variables Xn^m siid yn,m can be 
written in term of just one index, i.e. Xn and t/m, and we will just concentrate on the 
solution of the difference equation. 

3.1. Jacobi last multiplier on a two dimensional lattice 

Let us write an equivalent discrete expression defined on 4 lattice points of the linear 
first order partial differential equation ([T]) . We could write this expression in terms 
of the shift operators but we will use the difference operators (p7|) to follow closely 
the continuous case and to limit the number of points involved. Let us consider a two 
dimensional orthogonal lattice. In such a case ([T]) reads: 

As in the continuous case, if Un}m is a solution of fl3T]) . the 2x2 matrix 

A{Un,m,Un}m) _ ( ^xUn,ni AyUn,m 

has a determinant equal to zero if and only if the function Un,m. is also a solution of the 
difference equation (IM]) . The minors of the first row of (152]) are 

4(1) _ A ,,(1) 4(2) _ _ A ,,(1) 

and it is trivial to check by direct computation that A^^^ and A^'^^ satisfy the equation 

A^A^^l + AyA^^l = 0. (34) 

Then, as in the continuous case, there exists a function M^^rn, the Jacobi last multiplier, 
such that, 

4(1) _ A/f /■(!) 4(2) _ /M- f(2) (or) 

Consequently, the Jacobi last multiplier Mn,m satisfies a difference equation, which is 
the discrete analog of the differential one 0141) 

(A.M„,^tfi,„ + AyM^,^fi%^,) + AJi]l + Aji% = 0. (36) 



(32) 
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Given any particular, even trivial, solution M„ „i of the PAE fl36p . the solution of the 
overdetermined system fl35|) provides a solution of fl3T]) . We consider now a few examples 
of the calculus of the solution of linear difference equations using the JLM. 

3.2. Examples 

3.2.1. First example. We consider the equation: 

ym^x'^n,m ~l~ •^n^y'^n,m 0. (37) 

The equation for the Jacobi last multiplier ^ is the same as that for Un,m'- 

ym/^^M„,^rn + X„Aj^M„,^ = 0, (38) 

and a particular solution is obviously Mn,m = 1 which gives, taking into account ( l35l) . 
the following system of equations for Un^m- 

To solve this sytem of difference equations we need to specify the lattice. If we consider 
a uniform lattice in both variables i.e.: 

Xn+l - Xn = 5i, ym+l " Vm = ^2, (40) 

SO that Xn = Siu + xq and i/m = + yo, where xq and yo are arbitrary initial points, 
the system (139|1 reduces to a system of CAE's, one for each direction: 

Un+l,m = Un,m " - 6iXo, Un,m+1 = Un,m + + 52j/0- (41) 

Using the well know procedures for solving OAE's [8] we get a particular solution of 
f l37|) . depending on three arbitrary constants, i.e.: 

Un,m = M0,0 - ^{Xn + Xo){Xn - Xq - 6i) + ^{ym + ?/o)(2/m " 2/0 " (^2)- (42) 

3.2.2. Second example. We choose the equation: 

ymXn+l^xUn,m ~\~ •^Jij/m+l Aj^ttj^^m 0, ('^'^) 

which corresponds to (!3T]) with 

/n.m Vm^n+lj fn,rn ^nUra+l- (44) 

The equation for the Jacobi last multiplier can be written as: 



, Xn+2-Xn+l\+ ym+2 - ^n+l \ = . (45) 

A particular solution of the equation fllS]) can be obtained requiring that both its curly 
brackets be identically zero. In such a case we get a particular solution 

Mn,m = , (46) 

Vm+lXn+l 
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where a is an abitrary constant. If we introduce this resuh in ( l35l) with An}m and A< 
given by ( l33i) . we get the following system of compatible equations 



(2) 

n,m 



Ur. 



1 - 



Vrn+l 



(47) 



As in the previous example, we need the lattice equations to solve equations ( l47l) . 
Using again a uniform lattice in each variable, we get 



Xq + h^n, yrn = yo + hym, 



-ah^ I 1 



ahy I 1 



Xo + {n + 



[m + l)hy 



(48) 
(49) 



To solve the first equation we define, 
1 



ah. 



'^n,m ~t~ ri 



and the equation satisfied by Vn,m is: 

1 



Vn+l,' 



1 + ^ 



n 



(50) 



(51) 



which is the recursion equation for the Euler digamma function i/j. Then 

Xq 



(52) 



Substituting in the second equation in (149|) we easily obtain an equation fo Om- 



Solving this equation as in 
h, 



^11 

h^ 



2/0 



(m + l)hy 



hn, 



m 



(53) 



(54) 



where c is a constant. Then, the complete solution is 

Un,m = Uofi + a 



(55) 



4. Conclusions 



In this paper we have extended the results presented by Jacobi in 1844 to get solutions 
of linear partial differential equations to the case of partial difference equations in two 
independent variables. 
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If we consider an A^- dimensional lattice of independent coordinates, that is the 



lattice coordinates, a;'-*-', depend only on one index Ui, 



N, it is easy to see that 



we get into trouble as the minors are nonlinear functions and the difference operator, 
in contrast with the differential one, does not satisfy Leibniz rule. In fact denoting by 
Uni,...,nN the value of u at the point {xni , ■ ■ ■ , xi!^) and using the following notation: 

e, = (0,...,l(z),...,0), (56) 

(57) 



n = [Hi, . . .,nN), 
the discrete derivatives reads: 



[X. 



(1) 

rii ' 



,X. 



(JV)^ 



AiUr 



(i) 

rp ^ ' rp 



The difference equation is: 



N 



Ur 



0, 



(58) 



i=l 



where fn^ are some functions depending on a finite number of points in the lattice. If 
we know — 1 particular solutions of the equation, Un^ , ■ ■ ■ , u"n , we can construct 
the matrix 

^ AiUn ■ ■ ■ A^Un \ 



and define the minors An \ k — 
above (the column k is removed) 

^(fc) = (_i)'=-idet 



(59) 



y AiUn ■ ■ ■ L^N'an 

N, corresponding to the first line of the matrix 



k 



(60) 

\ Aiu''n~^^ ■ ■ ■ k ■ ■ ■ A^ui^"^^ / 

As in the continuous case (see Section [2]), the A^ — 1 solutions fix the coefficients of the 
difference equation (158|) up to a factor. Then, using Cramer's rule, we get 

AW = M^/W, k = l,...,N. (61) 

To find the compatibility relation we will closely follow the argument we used in 
the continuous case. Then, writing 

A^J:^ = det{Bu...,Bk,---,BN). (62) 

where Bi is the z-the column of fl60l) and 13^ means that the k-th column is absent, we 
get that the sum 

N 

J2 i-l)'~'det{B,,...,AkB„...,Bk,...,B^) (63) 

is equal to zero. In fact 

AkB, = (AfcA,n(,i), . . . , AfcA,n(f -1)) = A,S,, (64) 
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since A^Aj = AjA^ as the independent variables commute. Then 

det(i?i, . . . , AfcEj, . . . , Bk, . . . , B^) = det(-Bi, . . . , AjE^, . . . , Bf,, . . . , B^) 

= (-1)'=— 1 det(Si, . . . , . . . , A,5fc, ...,Bj^) (65) 

and the sum ( |63l) is zero. This is exactly the same equation we found in the continuous 
case. However, since Leibniz rule does not apply in the case of difference operators, the 
expression fl63l) is not equal to XlfeLi '^kAn\ as the difference operator does not follow 
the same rules as the differential operator when it is applied to a determinant. There 
are some additional terms, as A{fg) = fAg + gAf + h{Af){Ag), whose consequences 
have to be analyzed. They will be the content of a future work. 
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